home *** CD-ROM | disk | FTP | other *** search
/ Turnbull China Bikeride / Turnbull China Bikeride - Disc 1.iso / ARGONET / PD / MATHS / RLAB / RLAB125.ZIP / !RLaB / toolbox / tfd < prev    next >
Text File  |  1996-04-09  |  7KB  |  251 lines
  1. # function [y,min_t,max_t]=tfd(x,fa,sigma,threshold)
  2. #TFD Choi-Williams- and Wigner-Ville time-frequency distributions.
  3. #
  4. #  USAGE:    y=tfd(x,fa,sigma,threshold) 
  5. #            y=tfd(x,fa,sigma) 
  6. #            y=tfd(x,fa)       
  7. #            y=tfd(x)          
  8. #            [y,min_t,max_t]=tfd(...)
  9. #
  10. #            The default values for the non-defined parameters are
  11. #
  12. #            fa        = round(length(x)/2)
  13. #            sigma     = inf (i.e. a Wigner-Ville distribution)
  14. #            threshold = 0.001
  15. #
  16. #
  17. #  INPUT:
  18. #   x          Complex input signal (column- or row vector). If signal is
  19. #              real, use hilbert(x).
  20. #   fa         Length of frequency axis. (Default value: half the 
  21. #              signal length.)
  22. #   sigma      Degree of smoothing. The _smaller_ sigma, the _more_ 
  23. #              the crossterms are suppressed in the Choi-Williams 
  24. #              distribution. Nice value to start with: sigma=1.
  25. #              If sigma increases the result tends to a Wigner-Ville
  26. #              distribution. If sigma=inf (default) then an exact 
  27. #              Wigner-Ville distribution is returned.
  28. #   threshold  Minimum weightfactor (>0). In the Choi-Williams tfd the
  29. #              exponential weighting function in expands to
  30. #              infinity. However terms that are multiplied with a very low 
  31. #              weighting do add to the computation time, but do not add
  32. #              significantly to the final result. Using this parameter,
  33. #              only terms that have a weighting larger than "threshold" add
  34. #              to the result.
  35. #
  36. #  OUTPUT:
  37. #   y     -  Contains the distribution. Each row represents 
  38. #            a frequency, each column a time instant.  
  39. #   min_t -  first time-instant of distribution
  40. #   max_t -  final time-instant of distribution
  41. #
  42. #  
  43. # EXAMPLE 1: 
  44. # x=exp(2*pi*j* 6*((1:40)/40).^2);        # Make a chirp 
  45. # fa=20; sigma=inf;
  46. # [y min_t max_t]=tfd(x,fa,sigma);        # Make a WVD
  47. # contour((min_t:max_t),(1:fa), flipud(y))# Contour plot of WVD
  48. # hold on
  49. # plot((min_t:max_t),12*(min_t:max_t)/40) # Plot instantaneous frequency
  50. # hold off
  51. #
  52. # EXAMPLE 2:
  53. # figure(1)
  54. # x=exp(2*pi*j* 6*((1:80)/80).^2)+...     # Two 'crossing' chirps
  55. #   exp(2*pi*j*(1:80)/80.*(50-30*(1:80)/80)); 
  56. # fa=40; sigma=inf;
  57. # [y min_t max_t]=tfd(x,fa,sigma);        # Make a WVD
  58. # contour((min_t:max_t),(1:fa), flipud(y))# Contour plot of WVD
  59. # figure(2);sigma=1;
  60. # [y min_t max_t]=tfd(x,fa,sigma);        # Make a CWD
  61. # contour((min_t:max_t),(1:fa), flipud(y))# Contour plot of CWD
  62. # hold on
  63. # plot((min_t:max_t),12*(min_t:max_t)/80) # Plot instantaneous frequencies
  64. # plot((min_t:max_t),(50-60*(min_t:max_t)/80))
  65. # hold off
  66. #
  67. #  Author:   Rene van der Heiden, R.vanderHeiden@fel.tno.nl
  68. #  Ref:   -  H.I. Choi and W.J. Williams; IEEE Trans. Acoust., Speech, Sign.
  69. #            Proc., Vol. ASSP-37, pp 862-871, June 1989
  70. #
  71.  
  72. #  Translated to Rlab, 11/11/95, Ian Searle
  73.  
  74. #
  75. # This software may be freely used and modified for research and development
  76. # purposes. If you wish to use it for commercial gain please contact me. 
  77. # I provide ABSOLUTELY NO WARRANTY for this software.
  78. #
  79. #   Copyright: Rene van der Heiden
  80. #
  81.  
  82. tfd = function ( x, fa, sigma, threshold )
  83. {
  84.   local (x, fa, sigma, threshold)
  85.  
  86.   #
  87.   # Transpose signal if it is a column:
  88.   #
  89.  
  90.   NROW = x.nr; NCOL = x.nc;
  91.   if (NCOL == 1)
  92.   {
  93.     x=x.';
  94.     LENGTH=NROW;
  95.   else
  96.     LENGTH=NCOL;
  97.   }
  98.  
  99.   #
  100.   # Input check and, if necessary, assignment of default values
  101.   #
  102.  
  103.   def_sigma = inf();
  104.   def_fa = floor(LENGTH/2);
  105.   def_threshold = 0.001;
  106.   if (nargs == 1)
  107.   {
  108.     sigma = def_sigma;
  109.     fa = def_fa;
  110.     threshold = def_threshold;
  111.   else
  112.     if (nargs == 2)
  113.     {
  114.       sigma = def_sigma;
  115.       threshold = def_threshold;
  116.     else
  117.       if (nargs == 3)
  118.       {
  119.     threshold = def_threshold;
  120.       }
  121.     } 
  122.   }
  123.  
  124.   if (threshold <= 0)
  125.   {
  126.     disp("Parameter \"threshold\" is chosen less than or equal to zero. ");
  127.     disp("Using threshold="+num2str(def_threshold)+" instead.");
  128.     threshold = def_threshold;
  129.   }
  130.  
  131.   # In the next line we determine the size of the
  132.   # frequency window we will really use: window_use. (it should be uneven)
  133.   # 
  134.   # eg. if fa equals 8 then window_use becomes 7 and
  135.   # if fa equals 7 then window_use also becomes 7.
  136.   # 
  137.  
  138.   window_use = 2*floor((fa-1)/2)+1;
  139.  
  140.   #
  141.   # Calculate the number of timepoints at which the distribution will be 
  142.   # evaluated:
  143.   #
  144.  
  145.   length_cwd = LENGTH - window_use + 1;
  146.   
  147.   #
  148.   # Initialize cwd:
  149.   #
  150.  
  151.   cwd = zeros(length_cwd,fa);
  152.     
  153.   #
  154.   # Calculate the maximum size of the windows over which we smooth:
  155.   #
  156.  
  157.   wm = 2*floor((window_use-1)*sqrt(-log(threshold)/sigma));
  158.   
  159.   #
  160.   # Calculate for each tau and mu what the weights are:
  161.   #
  162.   
  163.   wgts = zeros((window_use-1)/2+1,wm+1);
  164.   mu = -wm/2:wm/2;
  165.   for (tau in 0:(window_use-1)/2)
  166.   {
  167.     if ((tau != 0) && (sigma != inf()))
  168.     {      
  169.       wgts[tau+1;1:wm+1] = exp(-mu.*mu*sigma/(4*tau*tau));
  170.     else
  171.       wgts[tau+1;1:wm+1] = (mu < 1);
  172.     } 
  173.   }
  174.  
  175.   #
  176.   # n is the time variable:
  177.   #
  178.  
  179.   min_t = (window_use+1)/2;
  180.   max_t = LENGTH - (window_use-1)/2;
  181.   for (n in min_t:max_t)
  182.   {
  183.   
  184.     #
  185.     # tau runs over the right part of window_use:
  186.     #
  187.  
  188.     for (tau in 0:(window_use-1)/2)
  189.     {
  190.     
  191.       #
  192.       # Summation over the second window to smooth the distribution:
  193.       #
  194.  
  195.       wm_tau = 2*floor(2*tau*sqrt(-log(threshold)/sigma));
  196.       mu_begin = -min(wm_tau/2,-1+    min(n+tau,n-tau));
  197.       mu_end   = min(wm_tau/2,LENGTH-max(n+tau,n-tau));
  198.       
  199.       weights = wgts[tau+1;mu_begin+wm/2+1:mu_end+wm/2+1];
  200.     
  201.       #
  202.       # Scale the weights:
  203.       #
  204.  
  205.       weights = weights/sum(weights);
  206.         
  207.       #
  208.       # And calculate the distribution:
  209.       #
  210.  
  211.       x1 = x[n+tau+mu_begin:n+tau+mu_end];
  212.       x2 = x[n-tau+mu_begin:n-tau+mu_end];
  213.       
  214.       cwd[n-(window_use+1)/2+1;tau+1] = sum(weights.*x1.*conj(x2));
  215.     }
  216.   }
  217.  
  218.   #
  219.   # Because the numbers in the CWD must be hermitic as a function
  220.   # of TAU, we can conjugate the found values of the CWD and copy
  221.   # them to the other half of the CWD.
  222.   #
  223.   
  224.   if (rem(fa,2) == 0)
  225.   {
  226.       
  227.     #
  228.     # if fa is even then write a zero in the cwd:
  229.     #
  230.          
  231.     cwd[n-(window_use+1)/2+1;fa/2+1] = 0;
  232.     for (tau in fa/2+2:fa)
  233.     {
  234.       cwd[;tau] = conj(cwd[;fa+2-tau]);
  235.     }
  236.     else
  237.     for (tau in (fa+1)/2+1:fa)
  238.     {
  239.       cwd[;tau] = conj(cwd[;fa+2-tau]);
  240.     }
  241.   }
  242.   
  243.   #
  244.   # And calculate the one dimensional simultaneous FT:
  245.   #
  246.  
  247.   y = real(fft(cwd',fa));
  248.  
  249.   return << y=y; min_t=min_t; max_t=max_t >>;
  250. };
  251.